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Abstract 

First principles electronic structure calculations based on density functional theory 
have been used to study the thermodynamic, structural and transport properties of 
solid solutions and liquid alloys of iron and oxygen at Earth's core conditions. Aims of 
the work are to determine the oxygen concentration needed to account for the inferred 
density in the outer core, to probe the stability of the liquid against phase separation, 
to interpret the bonding in the liquid, and to find out whether the viscosity differs 
significantly from that of pure liquid iron at the same conditions. It is shown that the 
required concentration of oxygen is in the region 25 — 30 mol percent, and evidence is 
presented for phase stability at these conditions. The Fe-0 bonding is partly ionic, but 
with a strong covalent component. The viscosity is lower than that of pure liquid iron at 
Earth's core conditions. It is shown that earlier first-principles calculations indicating 
very large enthalpies of formation of solid solutions may need reinterpretation, since 
the assumed crystal structures are not the most stable at the oxygen concentration of 
interest. 



1 Introduction 

The Earth's liquid outer core consists mainly of iron, but its density is about 10% too low 
to be pure iron (Birch, 1952), so that it must contain some light element. The nature of 
this element is still uncertain, and during the last 45 years the main candidates have been 
carbon (Birch, 1952; Clark, 1963; Urey, 1960; Wood, 1993), silicon (Birch, 1952; MacDonald 
and Knopoff, 1958; Ringwood, 1959, 1961, 1966), magnesium (Alder, 1966), sulphur (Clark, 
1963; Urey, 1960; Birch, 1964; Mason, 1966; Murthy and Hall, 1970; Lewis, 1973), oxygen 
(Dubrovskiy and Pan'kov, 1972; Bullen, 1973; Ringwood, 1977), and hydrogen (Birch, 1952; 
Fukai and Akimoto, 1983; Suzuki et al., 1989). For a given light element, it is also uncertain 
what concentration is needed to explain the inferred density in the core. The arguments for 
and against each of the candidate light elements have been reviewed by Poirier (1994). 

The aim of this paper is to use first-principles calculations to investigate the possibility 
that oxygen is the light element. First-principles calculations are well established as a reli- 
able way of predicting the thermodynamic, structural and dynamical properties of solid and 
liquid materials, including liquid metals (Stich et al., 1989; Kresse and Furthmuller, 1993). 
We have recently reported calculations of this kind on pure liquid iron under core conditions 



(Vocadlo et al., 1997; de Wijs et al., 1998), which show that it is a simple close-packed 
liquid with a viscosity not much greater than that of many liquid metals at ambient pres- 
sure, contrary to some earlier suggestions (Secco, 1995). We have also used first-principles 
simulations to investigate a liquid iron-sulphur alloy under the same conditions (Alfe and 
Gillan, 1998a). We showed that the properties of the liquid are scarcely affected by the small 
sulphur concentration needed to explain the observed density. 

The proposal that oxygen is the light element has a long and controversial history. Among 
the earliest proponents were Dubrovskiy and Pan'kov (1972) and Bullen (1973), the latter of 
whom suggested an outer core composition in the region of Fe 2 (equivalent to 12.5 wt%). 
Ringwood (1977) argued that oxygen should be seriously considered, and used seismic data 
to estimate the oxygen content as 28 mol percent (10 wt%). However, it is not completely 
certain whether the Fe/O liquid is thermodynamically stable against phase separation under 
Earth's core conditions at the composition that would be necessary to explain the density. It 
is known that the solubility of FeO in liquid Fe is very low (^ 1 mol percent) near the melting 
temperature of pure iron (1811 K) at atmospheric pressure (Distin et al., 1971). However, 
the solubility increases rapidly with temperature, becoming 6.5% at 2350 K (Fischer and 
Schumacher, 1978) and rising to the region of ~ 35 mol % at 2770 K (Ohtani and Ringwood, 
1984). According to Ohtani and Ringwood (1984), an extrapolation of the available phase 
measurements would suggest that the region of immiscibility disappears entirely above ~ 
3080 K at atmospheric pressure. It is also well established that the solubility of FeO in Fe 
increases with increasing pressure, and that the partial molar volume of FeO in liquid Fe 
is lower than that of pure liquid FeO itself. Ohtani et al. (1984) used their high pressure 
measurements on the solubility of FeO to suggest that the Fe/O system may show simple 
eutectic behaviour above a pressure of ~ 20 GPa, with no region of liquid immiscibility at any 
temperature. Subsequently, Ringwood and Hibberson (1990) showed by direct measurements 
that at 16 GPa addition of FeO to pure iron causes a depression of melting point, leading 
to a eutectic point at oxygen mole fraction of 28 % and a temperature of ca. 1940 K. 
Boehler's (1992) measurements on the melting of Fe/O mixtures are consistent with these 
ideas. Experiments of Knittle and Jeanloz (1991) and Goarant et al. (1992) on the reaction 
between lower mantle material and molten iron at pressures above 70 GPa revealed that 
the liquid dissolves significantly amount of FeO. However, Sherman (1995) has recently used 
first-principles calculations on crystalline Fe/O phases to argue strongly against significant 
amounts of oxygen in the core. His calculations gave values for the enthalpy of formation 
of crystals of composition Fe 3 and FeziO starting from Fe in the hexagonal close packed 
structure and FeO in the NiAs structure. The FeaO and Fe40 compositions were used 
to model substitutional and interstitial oxygen respectively. The enthalpies of formation 
were found to be so large that phase separation into FeO and Fe appears to be inevitable. 
However, it is not clear that Sherman's results have any relevance to the outer core, since 
in the liquid phase oxygen does not have to be either substitutional or interstitial. Even 
for the solid phase it is not obvious that Sherman's argument is robust, since Fe/O crystal 
structures other than those he studied might well give much lower enthalpies of formation. 

We are mainly concerned in this paper with the Fe/ O system in the liquid state. Our first- 
principles calculations, based on density functional theory and the pseudopotential method, 
will be used to address three questions: (a) Is the Fe/O liquid stable against phase separation 
under Earth's core conditions? (b) If it is, what oxygen concentration is needed to reproduce 
the observed density? (c) At this concentration, do the structural and dynamical properties 
of the liquid differ appreciably from those of pure liquid iron at the same pressure and tem- 
perature? Our first-principles simulations of the liquid will provide strong evidence that it is 
thermodynamically stable, and that the observed density requires an oxygen concentration 



of 25 — 30 mol-percent. In studying the properties of the liquid, we shall be particularly 
concerned with the viscosity, since this is one of the most poorly determined properties of 
the outer core, with estimates from different experimental and theoretical methods spanning 
many orders of magnitude (Secco, 1995). We shall also investigate a number of other prop- 
erties, including the nature of the short-range order, the atomic diffusion coefficients, and 
the electronic structure. 

Although we are mainly interested in the liquid, we shall also present some results for 
the energetics of various crystalline forms of the Fe/O system. These crystal calculations 
serve two purposes: first, they demonstrate that our techniques are in complete agreement 
with those used by Sherman in predicting large enthalpies of formation for FesO and Fe40 
in the structures he assumes; second, they demonstrate that there are crystal structures that 
give much lower formation enthalpies - an important fact in understanding how the Fe/O 
liquid can be stable against phase separation. It is not our intention to come to definite 
conclusions about the possible phase stability of the Fe/O solid solutions themselves, but 
our calculations suggest that their stability cannot be ruled out. 

The paper is organised as follows. In section 2, we summarise the first-principles tech- 
niques on which the work is based. Section 3 presents our calculations on the energetics of 
Fe/O crystals. In section 4 we report our results on liquid Fe/O, including its structural 
properties, the evidence for its stability against phase separation, and its dynamical and 
electronic properties. The final sections present discussion and conclusions. 

2 Methods 

In first principles calculations, the solid or liquid is represented as a collection of ions and 
electrons, and for any given set of ionic positions the aim is to determine the total energy and 
the force on every ion by solving the Schrodinger equation. This is a formidable task if the 
number of atoms is large, but it was made feasible by the introduction of density functional 
theory (DFT) many years ago (Hohenberg and Kohn, 1964; Kohn and Sham, 1965; Jones 
and Gunnarsson, 1989; Parr and Yang, 1989). DFT treats electronic exchange and corre- 
lation in a way that allows the electrons to be described by single-particle wavefunctions, 
with the interaction between them accounted for by an effective potential. DFT can be 
applied in two ways: all-electron calculations, or pseudopotential calculations. The first ap- 
proach includes such standard techniques as full-potential linearised augmented plane waves 
(FLAPW) and linearized muffin-tin orbitals (LMTO). In the pseudopotential approach, only 
valence electrons are explicitly treated, the effect of the core electrons being included by an 
effective interaction between the valence electrons and the cores. In both approaches, the 
accuracy with which the real material is described is governed by the approximation used for 
the electronic exchange-correlation energy. Until recently, the local density approximation 
(LDA) was the standard method. But in order to achieve the highest accuracy for transition 
metals it is essential to use an improved method known as the generalized gradient approxi- 
mation (GGA) (Wang and Perdew, 1991). The present work is based on the pseudopotential 
approach and the GGA. A non-technical review of first-principles calculations based on the 
pseudopotential approach has been given recently by one of the authors (Gillan 1997). 

There have already been extensive first-principles calculations on crystalline iron both 
at ambient pressure and at pressures going up to Earth's core values (Stixrude et al., 1994; 
Soderlind et al., 1996). The calculations have been performed using different all-electron 
techniques and the pseudopotential technique, and a variety of properties have been studied, 
including the equilibrium volume, the elastic constants, the magnetic moment, the volume as 



a function of pressure and lattice vibration frequencies. The agreement between results ob- 
tained with different techniques is generally very close, and the agreement with experimental 
data is also good. Particularly relevant here is the recent comparison of the pseudopotential 
results for the pressure-dependent volume of hexagonal-close-packed iron up to core pres- 
sures with earlier all-electron results and with experimental measurements (Vocadlo et al., 

1997) . 

Static first-principles calculations on crystals have been in routine use for many years. But 
to study liquids we need to do dynamical first-principles simulations, in which the calculated 
forces on the atoms are used to generate time evolution of the system, with every atom 
moving according to Newton's equation of motion. This kind of first-principles molecular 
dynamics (FPMD) pioneered by Car and Parrinello (1985) has been extensively used to 
study liquid metals (Stich et al., 1989; Kresse and Furthmuller, 1993; Holender and Gillan, 
1996; Kirchhoff et al. 1996a, 1996b), and it is known to give an accurate description of both 
structure and dynamics. 

The present work was performed mainly with the VASP code (Vienna Ab initio Simu- 
lation Package) (Kresse and Furthmuller, 1996a, 1996b). As usual in pseudopotential work, 
the electron orbitals are represented using a plane- wave basis set, which includes all plane 
waves up to a specific energy cut-off. The electron-ion interaction is described by ultra- 
soft Vanderbilt pseudopotentials (Vanderbilt, 1990), which allow one to use a much smaller 
plane-wave cut-off while maintaining high accuracy. When we perform FPMD with VASP, 
the integration of the classical equation of motion is done using the Verlet algorithm (1967), 
and the ground-state search is performed at each time-step using an efficient iterative ma- 
trix diagonalisation scheme and a Pulay mixer (1980). This method differs from the original 
Car-Parrinello technique, which treated the electronic degrees of freedom as fictitious dy- 
namical variables. In order to improve the efficiency of the dynamical simulation, the initial 
electronic charge density at each time step is extrapolated from the density at previous steps 
as described in our previous work (Alfe and Gillan, 1998a). In FPMD on metals, the dis- 
continuity of occupation numbers at the Fermi level can cause technical difficulties. Since 
we are interested in high temperatures in our liquid simulations, the electronic levels are 
occupied according to Fermi statistics corresponding to the temperature of the simulation. 
This prescription also avoids problems with level crossing during the ground state search. 
The FPMD simulations are performed at constant temperature (rather than at constant 
energy), using the Nose technique (1984). 

The iron pseudopotential we use is the same as that used in our earlier work (Vocadlo et 
al., 1997; de Wijs et al., 1998; Alfe and Gillan, 1998a,b), and was constructed using an Ar core 
and a 4s 1 3<i 7 atomic reference configuration. The oxygen pseudopotential was constructed 
using a He core and the 2s 2 2p 4 reference configuration. At Earth's core pressures the distance 
between the atoms becomes so small that the Fe(3p) orbitals respond significantly. The net 
effect is a small repulsion, which we determined from calculations on the h.c.p. crystal 
using a Ne core for the iron pseudopotential instead of Ar. We represent this repulsion by 
Fe-Fe and Fe-0 pair potentials, as in our earlier work (Vocadlo et al., 1997; de Wijs et al., 

1998) . The resulting corrections to energy and forces are generally small. Non-linear core 
corrections (Louie et al., 1982) are included throughout the work. 

For all the calculations to be reported, we use a plane-wave cut-off of 400 eV, which 
gives total energies converged within ss 10 — 20 meV/atom. In the calculations on crystals 
Brillouin-zone sampling is an important issue, and the sampling density we use is described 
in the next section. But for the liquid we use T point sampling, which experience suggests 
should be satisfactory. (We have done separate tests on pure liquid iron using 4 k-points, and 
we have found no detectable structural effects, while the average total energy difference with 



respect to the T point only calculations is of the order of 10 meV/atom, which is completely 
negligible for the purposes of the present work). The time step used in the dynamical 
simulations was 1 fs and we generally used a self-consistency threshold of 1.5 x 1CT 7 eV/atom. 
With these prescriptions the drift of the Nose constant of motion was less than xs 60 K per 
ps. 

3 Fe/O solid solutions 

We have used the techniques described in the previous section to calculate the equilibrium 
properties and the enthalpy of formation AH of various members of the Fe/O system, 
including those studied by Sherman (1995). Although Sherman's calculations and ours are 
both based on DFT, the technical methods are completely different, since Sherman used the 
FLAPW method, whereas ours are based on the pseudopotential approach. One of our main 
aims is therefore to make detailed comparisons with his results in order to ensure that the 
two methods agree about the energetics of the systems. We shall also point out that there 
are Fe 3 structures with enthalpies of formation much lower than those reported previously. 

We begin by presenting our results for the equilibrium density po (i- e - the density for 
which the pressure is zero), the bulk modulus K at this density, and the pressure derivative 
K' = dKj dP for crystals of Fe, FeO, Fe 3 and Fe40. Pure Fe is in the e structure (hexagonal 
close packed); FeO is in the B8 structure (the NiAs structure); FesO is in the structure 
obtained from face-centred cubic Fe by replacing the atoms at the corners of the conventional 
cube by O atoms; and Fe^iO is in the structure obtained from f.c.c. Fe by inserting an O 
atom at the centre of the conventional cube (see Fig. 1). We have done the calculations both 
spin restricted (the occupation numbers of every electronic orbital are equal for up and down 
spins, so that there are no magnetic moments) and spin unrestricted (the occupation numbers 
for up and down spins are allowed to vary independently). To sample the Brillouin zone (BZ) 
we have used Monkhorst and Pack grids (1976) using the sampling level that corresponds to 
20 and 36 k-points in the irreducible wedge of the BZ respectively for the cubic (FesO and 
Fe40) and the hexagonal (Fe(e) and FeO(B8)) crystals. Using these values, the total energies 
are converged within 5 meV/atom for the FesO, FeziO and FeO(B8) structures, and within 
10 — 15 meV/atom for the Fe(e) structure. We find that, except for Fe, all the structures 
are magnetic at low pressures, so that the system is stabilised if spin restriction is removed. 
This disagrees with the results of Sherman (1995), where a significant magnetic moment was 
found only for the FeO(B8) structure (a weak moment was found for FesO). 

Since we regarded the disagreement concerning magnetic properties as disturbing, we re- 
peated our calculations of the equilibrium magnetic moment of FesO and Fe^iO using a com- 
pletely independent electronic structure technique, namely the LMTO method (linearised 
muffin-tin orbitals), using the LMTO-46 code due to Krier et al. (1994). These calculations 
completely confirm the magnetic ordering in FesO and Fe40 and give numerical values for 
the magnetic moments that agree well with those given by our pseudopotential calculations. 
This suggests that the minimisation of the total energy with respect to magnetic moments 
may not have been fully under control in Sherman's work. 

For each system we have calculated the static internal energy E for a series of volumes, 
and fitted the results to the Birch-Murnaghan equation of state: 




where K is the zero pressure bulk modulus, K' = (dK/dP)p =Q , E is the equilibrium energy 
and Vq the equilibrium volume. Our calculated values of p , K and K' are compared with 
Sherman's results in Table 1, which also reports our calculated magnetic moments. The 
overall conclusion from the comparison is that the results agree well in the cases where 
magnetism is absent: pure Fe and spin-restricted Fe 3 and Fe40 (Sherman's calculations are 
effectively spin-restricted for Fe 3 and Fe40, since he found no moments). The equilibrium 
density agrees in those cases to better than 2%, and the values of K and K' to about 10%. 
However, our results show that magnetism has a strong effect on the equilibrium properties, 
so that it is important to treat it correctly. The case of FeO(B8) is problematic. Our spin- 
unrestricted calculations give a po value in respectable agreement with that of Sherman, 
but the agreement is poor for K and K' . It is clear that the accurate treatment of the 
volume dependence of magnetic moment is important in obtaining reliable values for these 
parameters, and our suspicion is that problems with moments may have affected Sherman's 
values. However, we shall stress below that magnetic effects become unimportant at core 
pressures, so the most important feature of Table 1 is the good agreement between the two 
sets of calculations for the non-magnetic cases. 

We turn now to the enthalpies of formation AH of Fe 3 and Fe 4 0, defined by: 

A#(Fe 3 0) = #(Fe 3 0) - #(FeO(B8)) - 2#(Fe(e)), 

A#(Fe 4 0) = H(Fe 4 0) - #(FeO(B8)) - 3#(Fe(e)), (2) 

where H = E + PV is the enthalpy per formula unit of each material. The total energy is 
taken directly from our Birch-Murnaghan fit, and the pressure P = —dE/dV is calculated 
from the derivative of the fitted form. We have done these calculations both spin restricted 
and spin-unrestricted, and our results are reported in the two panels of Fig. 2. In the spin- 
restricted panel we also show Sherman's enthalpy results. The most important conclusion is 
that AH for both Fe 3 and Fe 4 becomes very large at core pressures. In the range from 
135 GPa (core- mantle boundary) to 330 GPa (inner-core boundary) AH is at least 3 eV, 
which corresponds roughly to a temperature of 3.5 x 10 4 K, so that it is exceedingly unlikely 
that Fe 3 and Fe 4 could be thermodynamically stable in the assumed structures. 

Our spin unrestricted results indicate that the true values of AH are considerably lower 
than Sherman's results at low pressure, but at high pressure the differences between the 
magnetic and non-magnetic values of AH become very small, so the conclusion is unaffected. 
The detailed agreement with Sherman's values of AH is only moderately good, but again 
this does not affect the conclusions about the very large size of AH. 

We now want to ask whether the assumed crystal structures for Fe 3 and Fe 4 are 
actually the most stable. A glance at Wyckoff's book Crystal Structures (1964) shows that 
compounds having the composition A 3 B crystallise in a bewildering variety of structures. 
We have picked some likely candidates and calculated their formation enthalpy. Most turn 
out to be unfavourable, with AH values at least as great as those already reported in Fig. 2. 
However, we have discovered one that has a much lower value. This is the Bil 3 structure, 
which has a rhombohedral unit cell containing two formula units. Putting Fe 3 into this 
structure and relaxing both the atomic positions and the shape of the cell, we end up with 
the triclinic structure shown in Fig. 1. To characterise the structure briefly, we note that 
at 300 GPa each oxygen is surrounded by 11 Fe neighbours at distances of between 1.76 
and 2.57 A, and two O atoms at distances of 2.02 and 2.34 A (by contrast, in the cubic 
Fe 3 structure at the same pressure, each oxygen has 12 Fe neighbours at 2.05 A, and the 
nearest oxygen neighbours are at 2.9 A). We have calculated the fully relaxed total energy 
of this structure at several volumes, and the resulting structural parameters are reported in 



Table 2. Spin unrestricted calculations show that the structure is weakly magnetic, but the 
moment and the energy stabilisation are so small that the effects can be ignored. 

Fitting of the energies to the Birch-Murnaghan equation of state yields the enthalpy of 
formation shown in Fig. 2. Remarkably AH is very much lower than for the previous 
structures, and it decreases with increasing pressure. At the pressure of the inner core 
boundary it is only just over 1 eV. Since we arrived at this distorted Bil 3 structure in 
a rather haphazard way, it is quite likely that there are other Fe 3 structures with even 
lower enthalpies. We cannot say at present whether there are Fe 3 structures that are 
stable against phase separation under Earth's core conditions, but it certainly does not look 
impossible. The low AH for the Bil 3 structure will be highly relevant to our study of phase 
stability in the liquid Fe/O system. 

4 The liquid 

4.1 Thermodynamics 

Our aim in choosing the thermodynamic parameters for our liquid simulations was to model 
liquid Fe/O near the thermodynamic state it would need to have to reproduce the known 
outer-core density at the inner core boundary (ICB). The temperature at this point is very 
uncertain, with estimates ranging from 4000 to 8000 K (Poirier, 1991). We took the value 
of 6000 K, which is intended to be a reasonable compromise. However, the density and the 
pressure are quite accurately known to be ~ 12000 kg m -3 and m 330 GPa. This density 
is about 10% lower than it would be if the core consisted of pure iron (Birch, 1952). The 
main problem in choosing thermodynamic parameters is that we do not know in advance 
the required oxygen concentration, so that a certain amount of trial and error is needed. 

We started from our previous 64-atom simulation for pure liquid iron (Vocadlo et al., 
1997), which had a mass density of 13300 kg m~ 3 and a calculated pressure of 358 ± 6 GPa. 
Our first move was to hold the volume of the system fixed and to transmute the appropriate 
number of iron atoms into oxygen atoms to produce the density of 12000 kg m~ 3 . This 
resulted in a large reduction of the pressure, and we therefore reduced the cell volume to 
restore the original pressure. Naturally, this increased the density, and we therefore converted 
more iron atoms into oxygen to regain the density of 12000 kg m~ 3 . By repeating this cycle 
many times, one could in principle achieve the required density and pressure. But the 
calculations are very demanding, since at each state point one has to equilibrate the system 
and run it for long enough to obtain adequate statistics for the pressure, so that in practice 
a compromise between accuracy and computational effort is needed. After several iterations, 
we ended up with a simulation box containing 43 iron atoms and 21 oxygen atoms, i.e. mole 
fractions of xp c ~ 0.67 and xq 0.33. The resulting mass density of 11600 kg m™ 3 and 
pressure of 342 ± 4 GPa are close to the known values at the ICB. 

Since the mass density of 11600 kg ur 3 is slightly below the known value at the ICB, it 
is likely that the concentration of x ~ 0.33 is an overestimate. We have therefore taken a 
second thermodynamic state with a lower concentration. In order to facilitate comparisons 
with our calculations on crystalline Fe 3 0, we chose the value xo = 0.25. This second 
simulation was performed on a system of 48 iron atoms and 16 oxygen atoms at the mass 
density of 12200 kg m~ 3 , and the resulting pressure was 366 ± 8 GPa. We shall refer to the 
two simulations in the following as the '33% simulation' and the '25% simulation'. 

From the thermodynamic results just mentioned, we can estimate the oxygen concentra- 
tion that would be needed to reproduce the known density and pressure at the ICB. Interpo- 



lating between the calculated density values and applying a small correction for the slightly 
different pressures in the two simulations, we estimate that the mole fraction x<j = 0.28 
would reproduce the density 12000 kg m~ 3 at the ICB pressure. 

In the next sections we describe the structural, dynamical and electronic-structure prop- 
erties of the Fe/O liquid alloys. 

4.2 Structure 

We have simulated the 33% system for 4.2 ps after 2 ps of equilibration. The structural 
properties of the system have been inspected by looking at the partial radial distribution 
functions (rdf), g-F C F e (r), gF e o{r), and goo{ r )- The partial rdf's are defined so that, sitting on 
an atom of the species a, the probability of finding an atom of the species (3 in the spherical 
shell (r, r + dr) is 47rr 2 npg a p(r)dr , where rip is the number density of the species (3 (the mole 
fraction of species (3 times the total number of atoms per unit volume). 

We have calculated averages of the rdf's over different small time windows of the sim- 
ulation and we find no meaningful differences between the windows. This confirms that 
the system is well equilibrated. In Fig. 3 we display the rdf's calculated from the whole 
simulation. These show that the distance between neighbouring iron and oxygen atoms is 
significantly smaller than the iron-iron distance, the maximum of gFco( r ) being at ~ 1.7 A, 
while the maximum of gFeFc( r ) is at Pi 2.1 A. It is interesting to notice that goo( r ) has a 
first maximum at ~ 2.1 A, which is much greater than the chemical bond length expected 
for 0-0 single or double bonds (1.47 A and 1.21 A respectively). This is clear evidence 
that there is no covalent bonding between oxygen atoms. The presence of the 0-0 peak at 
rs 2.1 A indicates that oxygen atoms repel each other with an effective atomic diameter of 
w 2.1 A. This fact shows that oxygen has two effective sizes in the liquid: a small one when 
it interacts with iron and a large one when it interacts with itself. 

It is interesting to compare the structural properties of the alloy with those of pure liquid 
iron. In Fig. 4 we display the rdf calculated earlier for pure liquid iron at ICB conditions 
(Vocadlo et al., 1997) and the g FeFe calculated here. The two are not very different, the only 
apparent effect being the broadening of the peak in the liquid alloy, which is probably due 
to the greater disorder in the alloy. 

The integration of the first peak of the rdf's provides a definition of the coordination 
number N^p (the average number of neighbours of species (3 surrounding an atom of species 



where r c a a is the position of the minimum after the first peak of g a p. We find the values 
^FcFc = H O, -ZVpcO = 4.5, Nq Fb = 9.2, and Nq Q = 4.5. For comparison, the average 
coordination number found in our earlier simulation of pure liquid iron at ICB conditions was 
N§ eFe = 13.8 (Vocadlo et al., 1997). In interpreting these numbers, it is helpful to consider 
the coordination numbers that would be found if iron and oxygen atoms had exactly the 
same size and if atoms were packed in the same way as in pure liquid iron. In that case, the 
total number of neighbours of each iron atom, iVp eFe + iVp e0 , would be the same as in pure 
iron, whereas in fact it is 15.5. This increase of coordination number is clearly due to the 
smaller size of oxygen, which allows more atoms to be fitted into the first shell of neighbours. 
On the other hand, the total number of neighbours of each oxygen atom, Nq Fc + Nq Q is 
13.7, which is almost the same as the coordination number in pure iron. We interpret this 
as the result of two competing effects. The smaller size of oxygen would lead to a smaller 
coordination number if all atoms in its shell of neighbours were iron. But since on average 
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(3) 



4.5 of the neighbours are oxygen, which have a smaller size when interacting with iron atoms 
in the shell, the coordination number is increased again. 

We note that the structure of the liquid is very different from that of the cubic FeaO and 
Fe40 crystals discussed in section [| Oxygen atoms in these crystals have respectively 12 and 
6 iron neighbours. The coordination number of 9.2 in the liquid is roughly half way between 
the two. In the liquid, the radii from oxygen to iron and oxygen neighbours are equal to 
Ri 1.7 and ~ 2.1 A respectively, whereas in the cubic Fe 3 crystal at a similar pressure, 
the distances are 2.05 and 2.9 A. On the other hand, in the Bil 3 -structure FeaO, the O-Fe 
neighbour separation is spread over the range 1.76 — 2.57 A, and the 0-0 separation is in 
the range 2.02 - 2.34 A. 

We now want to ask whether our simulated system is really in a single phase and whether 
we can detect any sign of phase separation. In studying this it is very helpful to calculate 
the static structure factors S a p{k) defined by: 

Sap(k) = <p;(k) P/3 (k)>, (4) 

where (•) denotes the thermal average (in practice evaluated as a time average). Here, p a (k) 
is the Fourier component of the number density of species a at wavevector k, given by: 

p a (k) =iV Q - 1 / 2 5>xp(*k-r m ), (5) 
i=i 

where N a is the number of atoms of species a and r ai is the position of the ith atom of this 
species. Phase separation is associated with fluctuations of the concentrations of the two 
species, and the structure factors give us quantitative information about the intensities of 
these fluctuations. 

The connection between phase separation and structure factors can be made more precise. 
In the limit of zero wavevector, the structure factors of a liquid alloy can be rigorously 
expressed in terms of thermodynamic derivatives (Bhatia and Thornton, 1970): 



V,T. 



where \x a are the chemical potentials, and the notation indicates that the derivative is to 
be taken with the volume V, the temperature T and all chemical potential except \x$ held 
fixed. But the condition for thermodynamic stability with respect to phase separation is 

(<W^Mp,t > 0. (7) 

At the consolute point (the point in the phase diagram at which phases start to sepa- 
rate) the derivatives {dfi a /dx/3)p t T become zero. This implies that the matrix of derivatives 
(dfi a / 'dNp) v>T>N ' has vanishing eigenvalues, corresponding to variations of the numbers N a 
that maintain the pressure constant at fixed volume. But the matrix (dN a / djjLp) v,t,h' is the 
inverse of the matrix (dfi a /dNp)v,T,N' , so that when the latter becomes singular the former 
must acquire infinite eigenvalues. The consequence is that the values of the structure factors 
in the zero-wavevector limit must diverge if the system is unstable with respect to phase 
separation. This is also intuitively clear: as one passes from the miscible to the immiscible 
region, concentration fluctuations become ever larger, becoming of macroscopic size when 
the phases separate, and the increase in the fluctuations is reflected in the divergence of the 
quantities S a p{k — > 0). 



Our calculated structure factors for the 33% simulation are reported in Fig. 5. They have 
the form usually found in liquid alloys, with prominent peaks in SpeFe^) and Soo(k) in the 
region k ~ 4 A -1 signalling the approximate spatial periodicity associated with the packing 
of the atoms. The significant feature for present purposes is the lack of any anomalous 
behaviour at small wavevectors. We recognize, of course, that because of the limited size of 
the repeating simulation cell, there is a lower limit to the wavevector that we can examine, 
which in the present case is 0.86 A -1 . But at least in the accessible region of wavevectors 
there is no indication of any tendency towards phase separation. 

Before leaving the description of structure, we outline another method we have used to 
search for signs of phase separation. To explain this, let us imagine for a moment that the 
system had separated into phases of pure Fe and pure FeO. Then the Fe atoms in the Fe phase 
would have no oxygen neighbours, whereas the Fe atoms in the FeO phase might be expected 
to have 6 oxygen neighbours (we assume the FeO liquid to have a structure resembling that 
of crystalline FeO in the NiAs structure). On the other hand, in the unseparated Fe/O 
phase, the number of oxygen neighbours surrounding each Fe atom fluctuates around the 
value 4 — 5 (see above). We can therefore distinguish between the two situations by studying 
the probability distribution for the number of oxygen neighbours surrounding Fe atoms. To 
do this, we use the cut-off distance r° a p defined above to decide when an atom of species (3 
counts as a neighbour of an atom of species a, and we define the function P a p(n, as the 
probability that an atom of species a has n neighbours of species (3. If there is a complete 
phase separation, we expect PfcO to have peaks in the region of n = and n — 8, but if 
there is no separation we expect a single peak in the region of n = 4 — 5. Note that the 
rdfs contain less information than the P a p functions, and cannot by themselves deliver the 
discrimination we need. 

We present in Fig. 6 the function P Fe0 (n, Tp e0 ) calculated from our 33% simulation with 
the cut-off rp e0 = 2.5 A. This shows a single peak at n — 4, and no sign of any structure at 
lower or higher values of n. This means that there is no indication whatever of any tendency 
towards phase separation. 

4.3 Confirmation of phase stability 

Our failure to find any evidence of phase separation strongly suggests that the Fe/O liquid 
is in fact stable. But it might be objected that a simulation lasting only 4 — 6 ps does not 
give enough time for separation to occur, and that it would occur if the simulation were 
longer. To eliminate this possibility we have devised a method in which phase separation 
is artificially induced by an external force. We shall show that when the force is removed 
the phases spontaneously re-mix very rapidly. We performed this procedure on the 25% 
simulation (it is not significant that the 25% case was chosen for this, and we believe that 
the 33% system would have behaved in the same way). 

In our procedure, we notionally divide our cubic cell into two parts, consisting of the left 
region < x < 0.4 and the right region 0.4 < x < 1.0 (x is the coordinate along one of the 
edge directions in units of the cell length). The first step is to sweep all the oxygen atoms 
into the left region with an external force, so that this region contains something resembling 
FeO, while the right region contains pure Fe. To achieve this, we apply a constant force 
along the x-axis to all oxygen atoms lying in the right region. This force is in the positive 
x direction for 0.7 < x < 1.0 and in the negative direction for 0.4 < x < 0.7. No force is 
applied to the iron atoms, and these are left free to redistribute themselves. Initially, the 
magnitude of the force was taken to be 1 eV/A, but since this proved to be too weak it 
was increased to 3 eV/A. After ps 1 ps a complete phase separation was achieved, with all 



oxygen atoms in the left region, and we let the system equilibrate with the external force 
still present for a further 1 ps. We show in Fig. 7 a snapshot of one configuration taken 
from this period, which clearly shows the complete separation of phases. The external force 
was then switched off and the system was allowed to evolve for a further 2 ps. Remarkably, 
we found that after only 1 ps had elapsed the oxygen atoms became completely randomized 
throughout the cell, as can be seen in the snapshot shown in in the lower part of Fig. 7. 

To characterise these events quantitatively, we use the probability function P Fe0 (n, Tp e0 ) 
described in section |4.2| . Fig. 8 shows this function calculated by averaging over three short 
windows of 0.1 ps each, starting ps , 0.5 ps and 1.0 ps after the external force was switched 
off. The first window shows a clear bimodal form, as would be expected for a two phase 
system. (We note that because our system is rather small, many atoms are near the boundary 
between the phases, so that the peaks in Ppeo are not as sharp as they would probably be 
in a larger system.) After 0.5 ps, the distribution has already become unimodal, and after 
1 ps it is very similar to what we showed in Fig. 6 for the equilibrium 33% simulation. The 
conclusion is clear: the system is not stable in the separated state and returns very quickly 
to the homogeneous state. 

In section |]we have used the calculations on Fe/O crystals to show that phase stability in 
high-temperature Fe/O systems might well be expected. In particular, we gave an example 
of an FeaO solid structure whose enthalpy of formation is only « 1 eV. We have used our 25% 
simulation to create another such structure. This was done simply by quenching the liquid 
at the rate of 3000 K ps" 1 until the atoms came to mechanical equilibrium in an amorphous 
structure (this is clearly a local minimum of the total energy function). We can regard this 
as a crystal with the FesO composition having an unusually large supercell. The calculated 
enthalpy of formation of this solid at the pressure of 290 GPa is reported in the left panel of 
Fig. 2, and we see that its stability is even greater that that of the Bil 3 structure proposed 
for FesO in section |3|. This confirms the idea that even more stable structures may yet be 
found. 



4.4 Dynamics 

In studying the dynamical properties of the Fe/O liquid, our main concern is with the 
viscosity. However, we first give results for the atomic diffusion coefficients D a , which give a 
simple way of characterising the motion of the atoms. These are straightforwardly calculated 
from the mean square displacement of the atoms through the Einstein relation (Allen and 
Tildesley, 1987): 

— (J2 l r m(*o + t) - r ai (t )\ 2 ) -»• 6D Q t, as t -> oo, (8) 

where r ia (t) is the position of the ith atom of species a at time t, N a has its usual meaning, 
and (•) is the thermal average, in practice evaluated by averaging over time origins to- In 
studying the long time behaviour of the mean square displacement, it is convenient to define 
a time dependent diffusion coefficient D a (t): 

I N a 

D a (t) = tttHE Ki(t + t)- r m (t )| 2 ), (9) 
otN a ^ 

which has the property that 

lim DJt) = D, 

t— >oo 



(10) 



In Fig. 9 we display the iron and the oxygen diffusion coefficients calculated using Eq. @. 
From this data we estimate D Fc m 0.8 x 10" 8 m 2 s" 1 and D Q m 1.0 x 10" 8 m 2 s _1 . These 
values should be compared with those obtained for pure liquid iron, D Fe « 0.4 — 0.5 x 10~ 8 
m 2 s -1 (Vocadlo et al., 1997; de Wijs, 1998), and those obtained for Fe and S in our Fe/S 
simulation D Fe « 0.4-0.6 x 10" 8 m 2 s" 1 , and D s « 0.4-0.6 x 10" 8 m 2 s" 1 (Alfe and Gillan, 
1998a). This means that the two species of atoms in liquid Fe/O diffuse somewhat more 
rapidly than the atoms in both liquid iron and the liquid Fe/S alloy at the same pressure 
and temperature. 

In the past, atomic diffusion coefficients have often been used to estimate the viscosity 
of liquids via the Stokes-Einstein relation, and this was the procedure used in our previ- 
ous work on liquid Fe and Fe/S. Since the diffusion coefficients are larger in the present 
case, this procedure would lead us to expect a lower viscosity for liquid Fe/O. We have 
recently demonstrated (Alfe and Gillan, 1998b) that the viscosity can be more directly (and 
more rigorously) calculated in first-principles simulations using the Green-Kubo relations, 
i.e. the relations between transport coefficients and correlation functions involving fluxes of 
conserved quantities (Allen and Tildesley, 1987). The shear viscosity rj is given by: 



where V is the volume of the system and P xy is the off-diagonal component of the stress tensor 
P a /3 (a and (3 are Cartesian components). The stress tensor is straightforwardly calculated, 
so that the stress autocorrelation function (SACF) (P xy (to + t)P xy {to)) can also be obtained, 
but at first sight it might appear that very long simulations would be needed to gather 
adequate statistical sampling. However, we have recently shown that perfectly adequate rj 
values can be obtained with surprisingly short runs, and we have reported results for liquid 
aluminium and liquid Fe/S (Alfe and Gillan, 1998b). 

In the left panel of Fig. 10 we display the the average of the five independent components 
of the traceless SACF divided by its value at t — which we denote by <f>(t). Since the 
traceless part of the stress tensor has zero average, <p(t) goes to zero as t — > oo. The 
statistical error on <f>(t) for all values of t is ~ 5% of the value at t = 0, and after 0.2 — 0.3 
ps the magnitude of (f)(t) falls below that error. In the right panel of Fig. 10 we display 
the integral J * dt'(j){t') of <f){t) as a function of time. The limiting value of the integral for 
t — > oo is the shear viscosity. The error that one makes in evaluating that integral grows 
with time, since one integrates the noise together with <p(t). We estimated the error in the 
integral as a function of time using the scatter of the SACF's. Combining this estimate 
with an analytic expression for the error, we obtain the error estimate displayed in Fig. 10. 
From the point where <f>(t) falls below the noise one integrates only the latter, so one gains 
nothing by evaluating the integral beyond that point. If we assume that 4>(t) is zero above 
t ~ 0.3 ps, we obtain the value rj = 4.5 ± 1.0 mPa s. This value is roughly half the value 
reported earlier for liquid Fe (Vocadlo et al., 1997) and Fe/S (Alfe and Gillan, 1998a) at ICB 
conditions, and is not very much greater than the viscosity of typical liquid metals under 
ambient conditions; for example, the viscosity of liquid aluminium at atmospheric pressure 
100 K above its melting point is 1.25 mPa s (Shimoji and Itami, 1986). 

This result confirms our earlier conclusion (de Wijs et al., 1998) that the viscosity of the 
outer core is towards the lower end of the wide range of theoretical and experimental values 
reported in the literature. 




(11) 



4.5 Electronic structure 



We have studied the electronic structure of our simulated Fe/O liquid in order to shed light 
on the nature of the bonding between the atoms and to help to interpret the structure 



discussed in section [4.2| . The main tools used in this analysis are the electronic density of 
states (DOS) and the local density of states (LDOS). The DOS represents the number of 
electronic states per unit energy as a function of energy, while the LDOS is a projection of 
the DOS onto states of chosen angular momentum on atoms of chosen species. In performing 
this projection we took spherical regions of radius R on the atoms, and in practice we chose 
R = 0.6 A for both iron and oxygen. This R value is considerably smaller than half the 
interatomic distance in the liquid (see Fig. 3), so we expect to distinguish clearly between 
the electronic structures on different atoms. The results are not averaged over time but are 
calculated from the electronic energies and wavefunctions for selected time steps taken from 
the 25% simulation. 

The calculated DOS is shown in the upper panel of Fig. 11, and consists of four main 
features (energies are referred to the Fermi energy): two fairly narrow peaks at —24 and 
— 11 eV; a large dominant peak spanning the range —9 eV to 3 eV; and a broad feature 
extending well above the Fermi energy. The LDOS shown in the lower panel of the Figure 
allows us to interpret these features. The peak at —24 eV consists entirely of 0(2s) states. 
The peak at —11 eV is mainly 0(2p), but with an appreciable contribution of Fe(3d). The 
dominant peak from —9 eV to 3 eV consists mainly of Fe(3cf) states, but there is a significant 
peak associated with 0(2p) states just above the Fermi energy. The broad feature above the 
Fermi energy comes from both Fe(3c?) and 0(2p) states. 

To see the implications of this structure for the interatomic bonding, we note that if 
the 0(2p) states were much lower in energy than the Fe(3<i) states, there would be very 
little hybridisation between the two kind of states, the 0(2p) levels would be completely 
filled, and the oxygen atoms would carry a net charge of — 2|e|. The Fe-0 bonding would 
be ionic. On the other hand, if the 0(2p) states were at the same energy as Fe(3d), we 
should expect strong hybridisation and little charge transfer, so that the bonding would be 
covalent. It is clear from Fig. 11 that the Fe-0 bonding is intermediate between ionic and 
covalent. The 0(2p) states are below Fe(3d), but not low enough to suppress hybridisation. 
There is a clear splitting of the 0(2p) levels into bonding and anti-bonding orbitals, though 
the bonding orbital clearly has much larger weight (in pure covalent bonding we should 
expect the weights to be equal). The peak in Fe(3c?) at —11 eV is also clear evidence for 
hybridisation between 0(2p) and Fe(3d). The implication is that there is a partial charge 
transfer from Fe to O, but not enough to give oxygen a charge of — 2|e|. 

The bonding between Fe atoms is metallic. The partial filling of the Fe(3<i) levels gives 
the well known bonding mechanism emphasised by Friedel's (1969) analysis of the cohesive 
and elastic properties of transition metal crystals. 

Since the oxygen atoms carry partial charges and their 2p orbitals are almost full, no 
covalent bond is expected between them, and this is also clear from Fig. 3. To investigate 
the electronic state of oxygen in more detail, we have calculated the LDOS for oxygen atoms 
in different environments. We show in Fig. 12 the LDOS for two oxygen atoms denoted O a 
and Ob, which have been chosen from the 25% simulation so that O a has 10 Fe neighbours 
and one O neighbour while Ob has seven Fe and four O neighbours. If there were any 
covalent bonding between O atoms, we should expect a larger bonding-antibonding splitting 
of the 2p states for the Ob atoms than for O a , and we might expect a similar splitting (or 
at least a broadening) of the 0(2s) states for Ob- The LDOS curves show neither of these 
effects. Instead, the main difference is the upward shift of the peaks for Ob compared with 



O a . We believe this is direct evidence for the partial charge transfer: the valence electrons 
of Ob feel a repulsive electrostatic potential due to the partial negative charges on the four 
oxygen neighbours, which raises their energy. 

The main bonding mechanisms in the liquid are therefore the ionic-covalent Fe-0 bond 
and the metallic Fe-Fe bond. Since the 0(2p) atomic orbitals are more compact than Fe(3d) 
orbitals we expect the Fe-0 distance to be shorter than the Fe-Fe distance, and this effect is 
clear from the rdf 's shown in Fig. 3. The partial charge transfer presumably also contributes 
to the shortening of this distance. 

5 Discussion 

Two of our main aims in this paper have been to probe the phase stability of liquid Fe/O 
under Earth's core conditions, and to determine the oxygen concentration that would be 
needed to reproduce the known core density. In practice, these aims must be taken to- 
gether: we want to know whether the alloy is thermodynamically stable at the appropriate 
concentration. 

The results we have presented leave little doubt that if the liquid is stable then the 
mole fraction of oxygen must be in the region 25 — 30% (our best estimate is 28%), because 
anything much less than this would give a pressure that is too low at the known density. 
This is essentially the same as the value proposed many years ago by Ringwood (1977). 
In judging the robustness of this conclusion we recall some important facts: First, DFT 
electronic-structure methods of the type used here generally predict the density of materials 
at a given pressure to within a few percent. Particularly relevant here are recent DFT 
calculations (Stixrude et al., 1994; Soderlind et al., 1996; Vocadlo et al., 1997) on h.c.p. iron 
over the pressure range — 350 GPa which are in excellent agreement with each other and 
with the experimental results. Similar comparisons for FeSi (Vocadlo et al., 1997) are also 
relevant. DFT calculations on oxides (including transition-metal oxides) generally predict 
the density with similar accuracy. A second important fact is that our earlier first principles 
simulations of pure liquid iron, based on exactly the same techniques, gave a prediction of 
the density at the pressure of the inner core boundary which is correct to ~ 2%. (In fact 
the comparison was done the other way round: at the density of 13300 kg m -3 and the 
temperature of 6000 K, our simulations gave a pressure of 358 GPa, compared with the 
value of 330 GPa estimated from experimental data.) Third, our estimate of the oxygen 
concentration is based on changes of pressure and density compared with our simulated pure 
iron. The calculations should be even more reliable for these differences then they are for the 
absolute values. We therefore believe that our value for the oxygen concentration required 
should be subject to an error of no more that as 5%. 

It is interesting to compare the liquid composition with that for the case of Fe/S. In 
our recent simulations (Alfe and Gillan, 1998a), we showed that liquid Fe/S is not far from 
reproducing the known pressure and density at the inner core boundary with a sulphur mole 
fraction of 18%. At this composition, the mass density of 12330 kg m~ 3 gave a calculated 
pressure of 349 ± 6 GPa. If we make a correction to bring the density to 12000 kg m~ 3 , we 
find a sulphur mole fraction of 23%. This means that to achieve the required reduction in 
density we actually need a higher mole fraction of oxygen than of sulphur, even though the 
atomic mass of oxygen is only half that of sulphur. The reason is that the oxygen atom is 
smaller, a point to which we return below. 

The question of phase stability is more complex. What seems certain is that earlier 
arguments against phase stability based on ab initio calculations of the energetics of FesO 



and Fe40 crystals (Sherman, 1995) do not really deliver the intended conclusion. This is not 
because the calculations were wrong. On the contrary, our calculations fully support their 
correctness. It is simply that the crystal structures assumed were not the most stable. We 
have presented an alternative structure for FesO which gives a formation enthalpy that is 
low enough 1 eV) to make phase stability under core conditions quite plausible. 

We have used our simulations to probe the phase stability of the Fe/O liquid in the appro- 
priate region of concentration, and all the indications are that it is stable. We therefore fully 
support the conclusion that has been drawn from high pressure experimental measurements 
(Ohtani et al., 1984; Ringwood and Hibberson, 1990) that liquid Fe and FeO are miscible 
under core conditions in the concentration region of interest. However, a word of caution is 
in order. The fact is that our simulated systems are rather small, and it is quite conceiv- 
able that a system that would be unstable in the thermodynamic limit could be stabilised 
by the artificial periodic boundary conditions in small simulation cells. Nevertheless, our 
calculations certainly provide strong support for the thermodynamic stability of liquid Fe/O 
under core conditions at the relevant concentration. In the end, we believe that the definitive 
theoretical approach to this question will be first principles free energy calculations on the 
appropriate solid and liquid phases. These would be computationally very demanding, but 
should certainly be feasible in the near future. (First principles free energy calculations have 
recently been used with success to calculate the melting properties of silicon and aluminium 
(Sugino and Car, 1995; de Wijs et al., 1998).) 

The small size of oxygen is clear from our analysis of the liquid structure: the Fe-0 
nearest neighbour distance (the position of the first peak in the rdf) is only at pa 1.7 A, 
compared with k 2.1 A for the Fe-Fe distance. We recall that in the Fe/S liquid, the Fe-S 
distance was pa 1.95 A (Alfe and Gillan, 1998a). An important feature of the liquid is that 
each oxygen has on average only 9 iron neighbours, whereas it has between 4 and 5 oxygen 
neighbours. This atomic environment of oxygen is very different from that produced by 
the cubic structure of FesO. By contrast in the Bil 3 structure of FesO oxygen has 11 iron 
neighbours and 2 oxygen neighbours. The dilemma in making crystal structures for Fe/O 
solid solutions at core pressures is that we want to achieve close packing because of the high 
pressure, but the atoms that are being packed have different sizes. It seems that the cubic 
structure of FesO is not a good solution. The B1I3 structure is better, even though it means 
putting oxygen atoms next each other. There may be yet better ways. 

Our analysis of the electronic structure of the liquid gives further insight. Here, the 
important feature is that Fe-0 bonding is only partially ionic, with a substantial covalent 
contribution. The implication that there is partial electron transfer from Fe to O is relevant, 
because presumably if O carried a full ionic charge of — 2|e| oxygen atoms would be more 
reluctant to become neighbours of each other. The electronic structure shows that there is 
no detectable covalent interaction between oxygens, so the fact that they become neighbours 
cannot be attributed to covalency. 

Finally, we have studied the diffusion coefficient of iron and oxygen atoms and the viscos- 
ity of the liquid at 33% composition. The finding is that both species diffuse more rapidly 
than in either pure liquid iron or the liquid Fe/S alloy: the diffusion coefficients have roughly 
twice the value that they have in pure Fe and Fe/S at the same pressure and temperature. 
We also find that the viscosity of the Fe/O liquid is about half what it is in those other 
systems. This means that if the major light element in the outer core was oxygen, our 
earlier conclusion (de Wijs et al., 1998) about the low viscosity of the outer core would be 
confirmed, and indeed strengthened. 



6 Conclusions 



We are led to the following conclusions: if the light impurity in the outer core is mainly 
oxygen, then its molar concentration would have to be ~ 28%; in this region of concentration, 
we have strong evidence that the liquid is stable against phase separation; the proposed 
miscibility is not in conflict with the large formation enthalpies predicted by earlier ab initio 
calculations, because those calculations were based on assumed structures that are not the 
most stable; the proposed Fe/O liquid alloy has an even lower viscosity than that of pure Fe 
and the relevant Fe/S alloy under the same conditions. 
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Table 1: Calculated equilibrium properties of crystals in the Fe/O system: equilibrium 
mass density po (kg m~ 3 ), bulk modulus K (GPa), the pressure derivative K' = dK/dP, 
and magnetic moment per atom /j, (units of Bohr magneton). Results are given for the 
present spin unrestricted and restricted pseudopotential (PP) calculations and the FLAPW 
calculations of Sherman (1995). The structures of the cubic Fe 3 0, Fe 4 and Bil 3 -structure 
Fe 3 are shown in Fig. 1. 
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Table 2: Cell parameters and pressure P of Fe 3 in the Bil 3 structure (see Fig. 1) calculated 
at a series of volumes (per Fe 3 unit). The quantities a, b, c are the magnitudes of the 
primitive translation vectors, and a, (3, 7 are the angles between the pairs (a,b),(a,c) and 
(b, c) respectively. 



Figure captions 



Fig. 1: The crystal structures of cubic Fe 3 (left), cubic Fe40 (centre) and the Bil 3 form of 
Fe 3 (right) used to calculate the formation enthalpies of Fe/O solid solutions. Large and 
small spheres represent iron and oxygen respectively. 

Fig. 2: Calculated enthalpies of formation (per formula unit) of solid solutions having com- 
positions FesO and Fe40. Left panel shows spin-restricted results from present work com- 
pared with FLAPW results of Sherman (1995); right panel shows present spin-unrestricted 
results. Key to style of curves: present cubic Fe 3 , cubic Fe 3 of Sherman (1995) — , 

present cubic Fe 4 ■ - cubic Fe40 of Sherman (1995) , present Fe 3 in the Bil 3 

structure - ■ ■. The isolated point shows formation enthalpy of the amorphous structure 

obtained by quenching the liquid (see Sec. (|). 

Fig. 3: Radial distribution functions g a p{T) obtained from simulation of liquid Fe/O at 
oxygen molar concentration of 33%. 

Fig. 4: The iron-iron radial distribution function ^FeFe(^) from the present simulation of 
liquid Fe/O (oxygen molar concentration of 33%) compared with g-peFe{ r ) from simulation of 
pure liquid iron at similar pressure and temperature (Vocadlo et al., 1997). 

Fig. 5: Partial structure factors S a p{k) calculated from simulation of liquid Fe/O at oxygen 
molar concentration of 33%. 

Fig. 6: Probability distribution Pp e o(n, rp e0 ) for number n of oxygen neighbours of an iron 
atom calculated from simulation of liquid Fe/O at oxygen molar concentration of 33%. 

Fig. 7: Snapshots of the simulated liquid Fe/O system along three principal Cartesian axes. 
Top three panels show a configuration from the period when phase separation was artificially 
induced by application of an external force. Bottom three panels show a configuration from 
the later period after removal of the external force. Large and small spheres represent iron 
and oxygen respectively. 

Fig. 8: Probability distribution Pp e o(n, r FeQ ) for number n of oxygen neighbours of an 
iron atom calculated from simulation of liquid Fe/O at oxygen molar concentration of 25%. 
Results are average values for three windows of length 0.1 ps at times of ps, 0.5 ps and 1 
ps after removal of the external force used to induce phase separation. 

Fig. 9: Time dependent diffusion coefficients D a (t) for iron and oxygen calculated from the 
simulation of liquid Fe/O at oxygen molar concentration of 33%. 

Fig. 10: Average over the five independent components of the autocorrelation function 
of the traceless stress tensor <j>(t) (left panel) and viscosity integral (solid curve) with its 
statistical error (dashed curve) (right panel). 

Fig. 11: Electronic density of states (upper panel) and local densities of states (lower panel) 
calculated for liquid Fe/O at oxygen molar concentration of 25%. Energy is referred to the 
Fermi energy Ef. 

Fig. 12: Local densities of states for two selected oxygen atoms taken from the simulation 
of liquid Fe/O at oxygen molar concentration of 25%. Atom O a has 1 oxygen neighbour 
and 10 iron neighbours; atom Ob has 4 oxygen neighbours and 7 iron neighbours. Energy is 
referred to the Fermi energy Ef. 
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